Computationally Efficient Technique for 
Nonlinear Poisson-Boltzmann equation 



Sanjay Kumar Khattri 

Department of Mathematics, University of Bergen, Norway, 
sanj ay@mi . uib . no 
http : //www. mi .uib.no/ ~sanjay 



Abstract. Discretization of non- linear Poisson-Boltzmann Equation equa- 
tions results in a system of non-linear equations with symmetric Jaco- 
bian. The Newton algorithm is the most useful tool for solving non-linear 
equations. It consists of solving a series of linear system of equations (Ja- 
cobian system). In this article, we adaptively define the tolerance of the 
Jacobian systems. Numerical experiment shows that compared to the 
traditional method our approach can save a substantial amount of com- 
putational work. The presented algorithm can be easily incorporated in 
existing simulators. 



1 Introduction 

Lets consider the following non-linear elliptic problem 

— div (e gradp) + f(p, x, y) = b(x. y) in Q and p{x,y)~p D on dflu ■ 

(1) 

The above problem is the Poisson-Boltzmann equation arising in molecular bio- 
physics. See the References |2I7I9I1UI11I12) . Here, fl is a polyhedral domain in 
M 2 , the source function b is assumed to be in L 2 (Q) and the medium property 
e is uniformly positive. 

A Finite Volume discretization of the nonlinear elliptic equation results in a 
system of non-linear equations 

F(p) := Ax p h + A 2 (p h ) - b h = . (2) 

Here, F = [Fi(p), F 2 (p), • • • ,F„(p)] T , Ai is the discrete representation of the 
symmetric continuous operator — div(egrad) and A 2 is the discrete representa- 
tion of the non- linear operator f(p,x,y). 

A Newton-Krylov method for solving the non-linear equation is given by 
the Algorithm^ In the Quasi-Newton method (see Algorithmic, we are solving 
the Jacobian equation (J(p^.) Zip = — F(pt)) approximately. We are solving the 
system J(pfe) Apk = — F(p^) + with ||rfe|| is chosen adaptively. The quasi- 
Newton iteration is given by the Algorithmic In the Algorithms ^ an d El II ■ \\l 2 
denotes the discrete L2 norm and maxit er is the maximum allowed Newton's 
iterations. It is interesting to note the stopping criteria in the Algorithmic We 
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Algorithm 1: Newton-Krylov Algorithm 
Mesh the domain; 
Form the non- linear system : F(p); 
Set the iteration counter : k = ; 

while k < maxit er or ||Z\p||i 2 < tol or ||F(p)||x, 2 < tol do 

Solve the discrete system : J{pk) Ap = — F(p k ) with a fixed 

tolerance; 

Pk+i = Pk + Ap; 

k++: 
end 



Algorithm 2: Quasi-Newton-Krylov Algorithm 
Mesh the domain; 
Form the non- linear system : F(p); 
Set the iteration counter: k = 0; 

while k < maxji er or \\Ap\\l 2 tol or ||F(p) || /,., tol do 

Solve the discrete system : J(pk) Ap = — F(pk) with a tolerance 

1.0 x lQ- (k+1) ; 

Pk+i = Pk + Ap; 

k++; 
end 



are using three stopping criterion in the Algorithms. Apart from the maximum 
allowed iterations, L 2 norm of residual vector (||F(p)||i 2 ) and also L 2 norm of 
difference in scalar potential vector (||Z\p||i 2 ) are being used as stopping criterion 
for the Algorithms. Generally in the literature, maximum allowed iterations and 
the residual vector are used as stopping criteria [91101111 and references therein] . 
If the Jacobian is singular than the residual vector alone cannot provide a robust 
stopping criteria. 

2 Numerical Experiment 

Let us solve 10 in the domain fl = [—1,1] x [—1,1] with k = 1.0 |2I7I8I9| . fi is 
divided into four equal sub-domains (see Figure ^) based on e. 

— V • (e Vp) + k sinh(p) = f in J? and p(x, y) = x 3 + y 3 on dQo ■ 

(3) 

For solving the linear systems, we are using ILU-preconditioned the Conjugate- 
Gradient (CG) method. For the Newton algorithm the tolerance of the CG 
method is 1.0 x 10~ 15 . For the quasi-Newton method the tolerance of the CG 
method varies with the iterations k of the Algorithm|21as follows : 1.0 x \0~( k+1 ) ; 
k = 0, 2, . . . , 14. Figures |21 El and |3 reports the outcome of our numerical work. 
The Figures |31 and 0] compares convergence of the quasi-Newton and Newton 
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methods. The Figure |21 reports computational complexity of the quasi-Newton 
and the Newton methods. It can be notice, even if initial iterations of the Newton- 
Krylov algorithm are solved approximately, the convergence rate of the algorithm 
remains unaffected. The Figure |21 shows that such an approximation saves a 
substantial amount of computational effort. 




3 Conclusions 

Quasi-Newton method for solving non-linear system of equation with symmetric 
Jacobian matrix is presented. Numerical work shows that the presented tech- 
nique is computationally efficient compared to the traditional Newton-Krylov 
method. An efficient solution technique for Poisson-Boltzmann equation is of 
interest to the researchers in computational chemistry, bio-physics and molec- 
ular dynamics. The presented algorithm can be easily implemented in existing 
simulators. 
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Fig. 3. Convergence of the Li norm of 
residual vector A(p). 
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Fig. 4. Convergence of the L2 norm of 
difference vector zip. 
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